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Abstract. We study the geometry of forces in some simple models for granular stack- 
ings. The information contained in geometry is complementary to that in the distribu- 
tion of forces in a single inter-particle contact, which is more widely studied. We present 
a method which focuses on the fractal nature of the force network and find good ev- 
idence of scale invariance. The method enables us to distinguish universality classes 
characterized by critical exponents. Our approach can be applied to force networks in 
other athermal jammed systems. 



1 Jammed Matter and Force Networks 

Aggregates of particles can be found in a disordered solid-like state resulting from 
the phenomenon of jamming [1-5]. Granular materials, colloidal suspensions and 
molecular liquids are but a few examples of such systems that present a non- 
zero yield stress while trapped in one of many accessible metastable states. If 
thermal fluctuations are irrelevant, the forces on each particle must balance. 
Each stable configuration is thus characterized by a highly irregular network of 
forces spanning the entire system. 

Experimental [6-10] and numerical [11-19] studies have identified two main 
distinctive features of these force networks. Firstly, strong fluctuations are found 
in the magnitudes of inter-particle forces. The associated distribution function 
P{F) displays two characteristic properties: (i) it decays exponentially at large 
forces and (ii) it exhibits a plateau or small peak at small forces, which has been 
identified as a signature of jamming. The second experimental observation is 
that large forces are concentrated along tenuous paths, which have been deemed 
"force chains" . While P{F) has been commonly used also as a characterization of 
these force chains, strictly speaking it provides no information about the spatial 
organization of forces. In fact, so far force chains have been identified mainly 
visually, and a quantitative characterization seems to be lacking. 

By drawing an analogy with percolation, in this Letter we develop a geomet- 
rical description which associates a set of critical exponents with an ensemble 
of force networks. We apply this approach to three different models of static 
granular media under uniform pressure. We find that they belong to different 
geometrical universality classes although P{F) displays similar features in all 
three of them. 
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2 Force clusters 

Consider an ensemble of configurations of a fixed number of jammed particles, 
obtained numerically or experimentally. Each configuration defines a contact 
graph G, where nodes correspond to particle centers and edges connect particles 
in contact. Assuming there is no friction, the inter-particle forces are normal 
to the particle surface, and the underlying force network can be represented 
by associating with each edge i of G the corresponding force magnitude Fi. 
To investigate the geometry of forces, rather then the underlying geometry of 
contacts, we choose a threshold / and look at the subgraph G{f) of G obtained 
by selecting only the edges with Fi > f. For / small, G{f) consists of a single 
connected component, but as / increases, G(/) breaks up into a number of 
disconnected clusters. An ensemble of force networks thus induces a family of 
probability distributions of cluster sizes p{.s, /) for different thresholds /, the 
cluster size s being defined as the number of edges in a cluster. 

If the forces were distributed independently for each i, e.g. uniformly be- 
tween and 1 , then the force clusters would simply be bond percolation clusters 
[20]. In that case, in the thermodynamic limit N — > oo, a phase transition oc- 
curs at a critical value fc of /: an infinite cluster exists with probability 1 for 
f < fc, and with probability for / > fc. At fc, the cluster sizes arc power-law 
distributed, p{s, fc) oc s""^, and the correlation length diverges as ^ oc |/ — fc\~'^ 
near the threshold. The scaling exponents r and u are universal, they are in- 
dependent of the underlying geometry, and in fact they do not depend on the 
local distribution of forces P{F) or even their correlations, as long as these are 
short-ranged. 

In an ensemble of force networks corresponding to a jammed system, force 
and torque balance on each particle cause dependence and long-range correla- 
tions between bonds. Nevertheless, if the average forces are uniform over the 
extent of the system, we expect to find a critical threshold fc and an associated 
set of universal scaling exponents. The analogy with percolation moreover sug- 
gests that these exponents are independent of P{F) and thus provide a new, 
complementary characterization of force networks. 

3 Criticality and finite size scaling. 

An efficient method is necessary to study the existence of scale invariance around 
the critical threshold from numerical and experimental data. While, strictly 
speaking, the system becomes scale-invariant only in the thermodynamic limit, 
fc and the associated critical exponents can be extracted from data on systems 
of finite size using finite size scaling [21]. This describes the scaling of an ob- 
servable with the system size close to criticality: if a quantity X is expected to 
diverge as \ f — fc\~^ near fc in an infinite system, then in a system of size A'', 
it obeys the scaling law 



X{N,f) = N't'X{{f-fc)N'/''n 



(1) 
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Fig. 1. Examples of force networks (the thickness of the lines is proportional to the force 
magnitude) and corresponding force clusters close to the critical threshold: packing of 
400 grains in Model A (top) and packing of 200 grains in Model B (bottom). 



with d the spatial dimension and (j) = xldv. The scahng function X depends on 
a single rescaled variable x = {f — fc)N^^'^'^, and for a; ^ 1 it behaves as x~^, 
while for a; ^ it remains finite. 

Using measurements of X in systems of finite sizes, the parameters (f>, v and 
/c can be obtained from (1) in two steps. Assuming that X{N , /) as function of / 
displays a maximum Xm{N), from (1) the maxima for different A'' all correspond 
to the same maximum of X, hence XmiN) oc N'^. Plotting the amplitudes of 
the maxima versus N, we get the exponent </>. The values of fc and v can then 
be obtained by determining the best data collapse in the region around the 
maximum. 



4 Models studied 



Combining the finite-size scaling method with Monte-Carlo simulations, we stud- 
ied force-cluster criticality in three two-dimensional models of static granular 
matter under uniform pressure [22,23]. As all three models - which we will call 
A, B, and C for further reference - have been introduced earlier in other contexts, 
here we only define them briefly, without motivating in dc;tail their relevance to 
granular matter. In our view, they are are the simplest implementations of two 
fundamental ingredients of force networks, namely force balance on each grain 
and force randomness. 
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4.1 Snooker model 

To start with, we consider the "snooker-triangle packing" studied in [26,27]. It 
consists of a hexagonal packing of frictionlcss spherical grains confined within a 
triangular domain, with the same confining pressure applied on all sides of the 
triangle. A force network on this packing consists of repulsive forces in vectorial 
balance on each grain and consistent with the applied pressure. These constraints 
however do not define a single configuration of forces, but a whole set. Follow- 
ing Edwards' prescription [28], all such force networks are taken to be equally 
likely, similarly to a micro-canonical ensemble. We sample this ensemble with a 
Metropolis algorithm, using the parametrization of force networks developed in 
Ref. [29]. In Fig. 1 we show an example of a force network in this model and the 
corresponding force clusters for a threshold / = 0.94. 

4.2 Independent q-model 

We next consider consider the scalar g-model [24] , one of the first models intro- 
duced to account for the fluctuations of forces and appearance of force chains in 
a granular packing. Here we consider the massless g-model on a periodic tilted 
square lattice, which can be interpreted as a packing of rectangular bricks [25] . 
A uniform pressure is applied on the top of the packing and on each site a brick 
supports a weight Wij . Each brick transfers vertical forces i^/*"'^ and Fr^-'^ to its 
bottom left and right neighbors respectively. Vertical force balance is automat- 
ically satisfied by considering F^-^-'^ and Fr^-'^ respectively as fractions qij and 
1 — Qij of Wij, and randomness in force transfer is implemented by taking the 
qij uniformly distributed between and 1, independently for each site. Fig. 1 
shows a force network in this model and the corresponding force clusters for a 
threshold / = 0.7 (for unit external pressure). 

4.3 Microcanonic q-model 

Our third model is a variation on the g-model. We consider the same packing as 
in Sec. 4.2, but now, following Edwards' prescription, all allowed force networks 
- consisting of sets of vertical forces {{Fi''-'\ Fr^-'^)} - are equally likely. As 
shown in Ref. [30], this is equivalent to having the qij distributed with the joint 
probability distribution Yl^j Wij . The aim is to examine the influence of the form 
of the probability distribution by comparing independent and microcanonic q- 
models, and the difference between scalar and vectorial conservation laws by 
comparison with the snooker model. 

5 Results 

A convenient observable to study is the second moment of the distribution of 
cluster sizes, (s^(A^, /)) = / s^p{s,fc), where the contribution from the largest 
cluster in each configuration is omitted, and the system size N is defined as the 
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total number of edges. In all three models defined above, we find that {s'^{N, /)) 
displays a maximum as function of / for fixed N. The amplitudes of the maxima 
as functions of N follow sharp power- laws shown in Fig. 2 (a), thus confirming 
the existence of a critical threshold in each model. The corresponding critical 
exponent (j) is related via the hyper-scaling relation [20] to t, the exponent of 
the cluster-size distribution at criticality, and D, the fractal dimension of the 
incipient cluster as = = D — 1. Higher moments {s"{N,f)) display a 
similar scaling with exponents (pn = , implying that the full distribution 

p{s, fc) approaches a scaling form around the critical threshold. 

The value of the critical threshold fc depends on the scale set by the external 
pressure. Under unit pressure, we found a different fc for each model. Fig. 3 
displays the scaling functions obtained by collapse of the data. The estimated 
values of the critical thresholds and exponents are summarized in Table 1, where 
the two-dimensional percolation exponents are also included for reference. 

In Fig. 2(b) we show the probability distributions P{F) of force magnitudes. 
In the independent g-model, P{F) is exactly exponential [31], while in the other 
two models it is exponential for large forces, and displays a peak at small forces. 




Fig. 2. Results of Monte-Carlo simulations for the three models defined in the text: (a) 
scaling of the maxima of the second moment {s'^{N,f)) of the distribution of cluster 
sizes (omitting the largest cluster in every configuration), as function of the total system 
size N; (b) probability distributions P{F) of force magnitudes, obtained from 100 
samples of systems of 10'* particles. 



Table 1. Values of the critical threshold fc and the critical exponents and ly obtained 
from Fig.2 and the data collapse shown in Fig. 3. For two-dimensional percolation, exact 
values are shown inside brackets. 

fc (t> = D~l u_ 

Independent g-model 0.7 ±0.01 0.69 ± 0.01 3.1 ± 0.1 

Snooker model 0.93 ± 0.01 0.89 ± 0.01 1.65 ± 0.1 

Microcanonic g-model 0.585 ± 0.05 0.81 ± 0.01 1.65 ± 0.1 

Percolation 0.895(43/48) 1.33(4/3) 
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Fig. 3. Scaling of the second moment (omitting the largest cluster in each configu- 
ration) as function of the threshold / and the system size N. (a)(s'^)(/) for different 
system sizes in the snooker model, (b-d) Data collapse obtained by expressing the 
rescaled second moment of cluster sizes N~'^ {s^) as function of the rescaled variable 
(/ — fc)N^^^'^ for the three models defined in the text: (b) snooker model, (c) indepen- 
dent g-model and (d) microcanonic g-model The values of the corresponding parameters 
fc, 4> a-nd v are summarized in Table 1. We do not show the data for very small system 
sizes where the collapse takes place only in a small region around the maximum. For 
the g-models, the systems studied had the same vertical and horizontal linear sizes. 



6 Discussion 



We have introduced a new approach to investigate the geometry of force net- 
works, based on statistics of clusters created by forces larger then a given thresh- 
old. The existence of a critical threshold uncovers a scale-invariance of force net- 
works, which we characterized by the critical exponents v and for the correla- 
tion length and the second moment of the cluster size distribution. In particular, 
in each network we identify a fractal object of dimension D = 1, given by the 
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incipient force cluster at the critical threshold. As shown in Table 1, we found 

three different sets of critical exponents for the three models we stiidicd, imply- 
ing that they belong to distinct geometrical universality classes, although their 
P{F) display similar features. Interestingly, for the snooker model, ^ is very 
close to the percolation value, but as the values of v are further apart and the 
scaling functions are different, it does not belong to the percolation universality 
class. 

Two distinct universality classes could have been expected a priori for the 
g-models on one hand and the snooker packing on the other. Indeed, the q- 
models are both directed and include only scalar conservation laws, while the 
snooker model is isotropic with vectorial conservation laws. The reason for the 
segregation of the independent and microcanonic in two different universality 
classes is more subtle. They differ only by the form of the probability distribution 
of forces, but in the independent case the distribution is Markovian from top 
towards bottom, while in the microcanonic case no such preferred direction of 
propagation exists. 

While in jammed matter the disorder in the underlying contact geometry 
plays an important role, we considered here only lattice models with fixed con- 
tact geometry. The force-cluster method can be applied in a straightforward 
fashion to ensembles of forces networks resulting from disordered contact net- 
works. By analogy with other critical phenomena, we however do not expect such 
randomness to modify the universality class. Indeed, the values we have found for 
u in combination with the Harris' criterion [32] suggest that geometric disorder 
is irrelevant. A further study indeed shows that introducing quenched disorder in 
the g-models does not modify the universality class, which in turn confirms that 
the scalc-invariance found here is in all aspects similar to equilibrium critical 
phenomena. 

While the results presented here clearly show that our method is able to 

discriminate between different scaling behaviors, a crucial question is whether 
any of the models belongs to the same universality class as a realistic two- 
dimensional system of grains under isotropic pressure. A recent study [33] of 
packings generated by molecular dynamics simulations showed that packings 
under isotropic pressure lead to the same scaling behaviour irrespectively of 
the applied pressure, the polydispersity of the grains, the coefficient of friction 
and the force law. Remarkably, the corresponding scaling exponents and scaling 
function appear to be the same as those obtained from the snooker packing. 

7 Outlook 

The existence of universality classes for force networks raises a number of new 
questions. First of all, what properties of a jammed system determine the uni- 
versality class of its force network? Our results suggest that that the isotropy 
of the applied force and the vector nature of the force balance are essential. On 
the other hand, packings under static shear might lead to another universality 
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class. Another relevant parameter could be the temperature in thermal systems 
which exhibit jamming, such as colloids. 

Our method based on force cluster criticality is clearly able to discriminate 
between the many models proposed for force networks [34,35,15,36-38]. In par- 
ticular, it shows that the Edwards' hypothesis, which proposes to consider all 
metastable states of a jammed system equally likely, leads quantitatively to the 
same scaling properties as found in force networks generated by molecular dy- 
namics simulations. 

Finally, the method developed here for force networks in jammed matter 
is clearly more general. It applies in principle to any ensemble of graphs with 
continuous variables on the edges, such as flux, transport or metabolic networks 
[39,40]. The corresponding universality classes could complement the topological 
characterizations of networks developed in the recent years [41] . 
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